//
//  CFDProblemFactory.cpp
//  NSESolver
//
//  Created by Patrick Zulian on 5/6/11.
//  Copyright 2011 ICS - USI - Lugano. All rights reserved.
//


#include "CFDProblemFactory.h"


namespace nsesolver {
    
    const Reference< CFDProblem > CFDProblemFactory::createProblem(const std::string &path,  const mpi::CartesianTopology & cart, const int imax, const int jmax)
    {
        using namespace std;
        
        const double alpha = 1;
        
        std::string pathNonConst(path);
        ExodusObjekt exodusObject(&pathNonConst[0]);
        
        
        Vector2d swCorner, neCorner;
        
        swCorner[0] = numeric_limits< double >::max(); 
        swCorner[1] =  numeric_limits< double >::max();
        
        neCorner[0] =  numeric_limits< double >::min(); 
        neCorner[1] =  numeric_limits< double >::min();
        
        // browse all coordinates in order to get maximal/minimal x and y values
        for (int i=0; i< exodusObject.getVarObj("coord")->range[0]; i++) {
            
            int nodeID = i;
            
            neCorner[0]  = max( neCorner[0], exodusObject.getVarObj("coord")->getDValueAt(nodeID*2));
            swCorner[0]  = min( swCorner[0], exodusObject.getVarObj("coord")->getDValueAt(nodeID*2));
            
            neCorner[1]  = max( neCorner[1], exodusObject.getVarObj("coord")->getDValueAt(nodeID*2+1));
            swCorner[1]  = min( swCorner[1], exodusObject.getVarObj("coord")->getDValueAt(nodeID*2+1));
               
        }
        
        
        NSE_SOLVER_ASSERT( neCorner[0] != swCorner[0] );
        NSE_SOLVER_ASSERT( neCorner[1] != swCorner[1] );
        
        
        double hX = (neCorner[0] - swCorner[0])/imax;
        double hY = (neCorner[1] - swCorner[1])/jmax;
        
        
        Obstacle obstacle( imax, jmax );
        
        assignObstaclesAndBoundary(obstacle, exodusObject, swCorner, neCorner, hX, hY);
        obstacle.assignFlags();
        
        BoundaryHandler bc(cart, obstacle);
        InitialValue iv(InitToOnes, InitToZeros);
        Reference< CFDProblem > refProblem = new CFDProblem( imax, jmax, alpha, swCorner, neCorner, bc, iv, cart );
        
        return refProblem;
    }
    
    
    void CFDProblemFactory::assignFlag(int el, int ssID, Obstacle & obstacle, ExodusObjekt & exodusObject, const Vector2d & sw, const Vector2d & ne, double hX, double hY ) {
        
      
        using std::min;
        using std::max;
        
        // This function receives an element and a sideset id
        // and assigns flags in the CFDProblemFactory::flag array employing the getIdx
        // routine.
        //
        // Here, maxX, minX, maxY, minY are the maximal dimensions of this element. For every quadrangle of
        // your discretization which is lying inside these dimensions will the flag changed to 
        // ssID. Based on this, you can set your boundary values but can also assign obstacles, as 
        // will be explained on the second exercise sheet.
        
        
        double maxX = -1e6, minX = 1e6, maxY = -1e6, minY = 1e6;
        ExVar *ss_prop = exodusObject.getVarObj("ss_prop1");
        ssID = ss_prop->getIValueAt(ssID-1);
        for (int i=0; i<4; i++) {
            
            int nodeID = exodusObject.getVarObj("connect1")->getIValueAt(el*4+i);
            
            maxX = max(maxX, exodusObject.getVarObj("coord")->getDValueAt(nodeID*2));
            minX = min(minX, exodusObject.getVarObj("coord")->getDValueAt(nodeID*2));
            
            maxY = max(maxY, exodusObject.getVarObj("coord")->getDValueAt(nodeID*2+1));
            minY = min(minY, exodusObject.getVarObj("coord")->getDValueAt(nodeID*2+1));
        }
        
        int iStart = (int) floor( (minX- sw[0])/hX);
        int iStop = (int)  ceil( (maxX- sw[0])/hX);
        
        int jStart = (int) floor( (minY-sw[1])/hY);
        int jStop = (int) ceil( (maxY-sw[1])/hY);
        
        
        for (int ii = iStart; ii<=iStop; ii++) {
            for (int jj = jStart; jj<=jStop; jj++) {
                
                double spY = sw[1] + (jj-0.5)*hY;
                double spX = sw[0] + (ii-0.5)*hX;
                
                if (minX < spX && spX < maxX && minY < spY && spY < maxY) {
                    if (ssID == 1) 
                        obstacle.addCell( ii-1, jj-1 );
                }
            }
        }
        
    }
    
    
    void CFDProblemFactory::assignObstaclesAndBoundary(Obstacle & obstacle, ExodusObjekt & exodusObject, const Vector2d & sw, const Vector2d & ne, double hX, double hY) {
        
        // if CFDProblemFactory stores an ExodusObjekt-Pointer, here for
        // each sideset (as assigned in cubit) flags will be asssigned.
        // This is basically done by calling assignFlag, which is called within the
        // sideset - loop of this function.
        
        
        
        ExVar* elementSideSetMember = exodusObject.getVarObj("elementSideSetMember");
        
        
        
        for (int el = 0; el< elementSideSetMember->range[0]; el++) {
            for (int si = 0; si< elementSideSetMember->range[1]; si++) {
                for (int ss = 0; ss< elementSideSetMember->range[2]; ss++) {
                    int ssID = elementSideSetMember->getIValueAt( el*elementSideSetMember->range[1]*elementSideSetMember->range[2] + si*elementSideSetMember->range[2] + ss);
                    if (ssID == -1)
                        continue;
                    
                    assignFlag(el, ssID, obstacle, exodusObject, sw, ne, hX, hY);
                }
            }
        }
    }
    
}
